Barkhausen instabilities from labyrinthine magnetic domains 
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Experimental investigations of the scaling behavior of Barkhausen avalanches in out-of-plane 
ferromagnetic films yield widely different results for the values of the critical exponents despite 
similar labyrinthine domain structures, suggesting that universality may not hold for this class of 
materials. Analyzing a phase field model for magnetic reversal, we show that avalanche scaling 
is bounded by characteristic lengthscales arising from the competition between dipolar forces and 
exchange interactions. We compare our results with the experiments and find a good qualitative 
and quantitative agreement, reconciling apparent contradictions. Finally, we make some prediction, 
amenable to experimental verification, on the dependence of the avalanches behavior from the film 
thickness and disorder. 



I. INTRODUCTION 

Understanding magnetic hysteresis is a fundamental 
theoretical problem with important practical implica- 
tions for magnetic devices^. Magnetic reversal is usu- 
ally associated with crackling noise — the Barkhausen 
effect — due to the interaction between magnetic domains 
and structural disorder^. The statistical properties of 
the Barkhausen noise are characterized by scaling laws, 
a signature of an undelying non-equilibrium critical be- 
havior, suggesting that it should be possible to separate 
magnetic materials into well defined universality classes 
characterized by the same values of the critical expo- 
nents and the same form of the scaling functions 4 . This 
program has been successfully carried over for bulk ma- 
terials with relatively simple parallel domain structures 
where it is possible to predict distinct universality classes 
depending on the sample microstructure (amorphous or 
polycrystalline/^. 

The situation is more complicated in thin films which 
can show complex domain structures due to the inter- 
play of magnetic anisotropics and dipolar interactions. 
In films with in-plane anisotropy, such as the MnAs films 
studied in Ref. |6j the domain structure crosses over 
from zigzag to rough as a function of the temperature 
and this change is reflected by a crossover in the uni- 
versality class of Barkhausen avalanche statistics^. The 
main effect of temperature in this material is to modify 
the strength of dipolar interactions^, a relevant param- 
eter for the critical behavior as confirmed by numerical 
simulations^. The magnetic properties of films with out- 
of-plane anisotropy have been intensively investigated in 
recent years due to their potentially higher bit packing 
density^. These systems typically display a labyrinthine 
domain structure as shown by high resolution magneto- 
optical techniques 9 or magnetic force microscop ^ 10 * 11 *. 
These techniques have been used to follow magnetic re- 



versal and measure Barkhausen avalanches, reporting ap- 
parently contradictory results for the critical exponent of 
the avalanche statistics. Whether universality holds for 
these systems remains thus an intriguing open question. 

Here, we use a phase field model describing magneti- 
zation dynamics in an out-of-plane film l 12 * 13 * to show that 
the scaling behavior of the avalanche size distribution is 
severely limited by dipolar interactions that set a char- 
acteristic lengthscale to the problem. The phase-field 
model overcomes the limitations of the dipolar random- 
field Ising model 14 which is plagued by lattice effects. At 
the same time the model, owing to its scalar structure, 
is much easier to simulate than the full micromagnetic 
equations. Our numerical results are in good agreement 
with experimental data and allow to clarify the origin 
of the observed discrepancies between the exponents re- 
ported in the literature. This paper is organized as fol- 
lows: in Sec II we illustrate the phase field model used 
for our calculations, in Sec. Ill we present the results of 
our calculation including hysteresis loops and avalanche 
statistics, in the last section we compare our results with 
available experimenatal data and we draw some conclu- 
sions. 



II. MODEL 

In the model, the magnetization M(R) at position R 
is perpendicular to the film plane and constant along its 
thickness d, as illustrated in the inset of Fig{l] There- 
fore to describe the domains behavior, we adopt a two- 
dimensional phase field model where a dimensionless 
scalar order parameter m(R, t) = M z (H,t)/M s repre- 
sents the local magnetization value on the film plane, 
M s being the saturation magnetization 12 . The energy 
per unit thickness is given by: 
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FIG. 1: (left) Hysteresis loops for different dipolar interaction strength 7. The inset sketches our model ferromagnetic thin film. 
Notice that the out-of-plane magnetization is constant along the film thickness. Panels (a-h) represent the film magnetization 
m(x,y) for different external field values along the hysteresis loop for 7 = 0.7 (a-d) and 7 = 0.5 (e-h). 
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the first term mimics the anisotropy energy with a dou- 
ble well potential V = -V (m 2 /2 - m 4 /4) givine; the 
order parameter two equivalent preferential orientations 
±1, with the constant Vo determining the barrier height, 
in particular Vo = K u /A with K u uniaxial anisotropy 
constant of the material. The second term represents the 
exchange interaction, opposing any magnetization vari- 
ation with strength A. The third term represents the 
long-range, non-local, dipolar interaction in the approxi- 
mation of small film thickness, given a local value of m, 
it promotes the magnetization reversal in its surround- 
ings (fio is the vacuum permeability). Pinning centers 
are modeled by a Gaussian uncorrelated random field 
H r with (ff r (R))=0 and (H r (R)H r (R')} = Atf(R-R'). 
Finally we consider a uniform external field H e that 
is slowly ramped up and down to produce an hystere- 
sis loop. We adopt a dimensionless scheme with unit 
length £ = ^ AjVo and unit field M s , defining r = H/£, 
h e = H e / M s , D = A/£ 2 M 2 and h r = H r /M s . The order 
parameter equation of motion can be derived in the ap- 
proximation of small time fluctuations, i.e. far from the 
critical point, through the hamiltonian functional deriva- 
tive dM(R,t)/dt = -A 6H[M(B.,t)]/SM(K,t) leading 
to: 
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& — Vo/fioM 2 and 7 = d/£47r are dimensionless pa- 
rameters and /jLqX is the unit time. For the numerical 
solution of eq. Q we used the finite- difference semi- 
implicit method of Ref. [12j All the simulations have 
been performed keeping a — 3 and the disorder strength 



D = 2 • 10 3 except when another value is explicitly spec- 
ified, the time step At = 0.5 and the squared mesh step 
Ax = Ay = 0.4. 



III. RESULTS 
A. Hysteresis loops and domain formation 

Typical hysteresis loops obtained by a slow ramp of 
the external field h e are shown in FigJT] together with 
selected plots of the domain structures (see also supple- 
mental movies 1 and 2). Depending on the strength of the 
dipolar term, tuned by the parameter 7, the magnetiza- 
tion dynamics can vary considerably. For relatively large 
7 (see Fig. [TJa-d)), we observe simultaneous nucleation 
of multiple droplet domains that grow as the external 
field is decreased below the saturation value. When the 
domain density is sufficiently high, domains merge giving 
rise to a labyrinthine structure. For smaller values of 7 
(see Fig. [TJe-h)) a single nucleated droplet grows form- 
ing a labyrinth. This process occurs suddenly at a critical 
value of the external field and corresponds to a jump in 
the hysteresis loop. Upon further decreases in the exter- 
nal field, domains expand by moving domain walls. It 
is remarkable that while the resulting labyrinthine do- 
main structures are similar the dynamics that generates 
them is completely different. The two types of dynam- 
ics have been observed experimentally in irradiated and 
non-irradiated Ni films at low temperature^. Fig{2] (a) 
and (b) describe the two different domain dynamics with 
a color scale highlighting subsequent avalanches. Both 
pictures have been taken in the central part of the hys- 
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FIG. 2: Panels (a) and (b) show the evolution of the mag- 
netization close to the center of the loop for 7 = 0.7 (a) and 
7 = 0.5 (b), the color scale highlights the shape of the do- 
mains at three different subsequent times. Panel (c) describes 
the bridging phenomenon showing the magnetization differ- 
ence for two consecutive times. Four avalanches are high- 
lighted having almost the same size, surrounded by smaller 
avalanches triggered by the bridging. Panel (d) shows the 
characteristic lengths in dimensionless units as a function of 
the ratio a/j. The dashed and dotted lines are the theoretical 
predictions: id = ot/j for the domain width and i w — y/2 for 
the domain wall width. 



teresis loop at the same values of the applied field for 
7 = 0.7 and 7 = 0.5 (see also supplemental movies 3 and 
4). In Figj2] (a), avalanches occur only close to the tails 
of the domain, leaving their width unchanged. When 
two neighboring domains are close to each other, a single 
avalanche bridges them together, yielding a characteristic 
avalanche lengthscale proportional to the domain width. 

In Figj2|b) avalanches of a wide variety of sizes oc- 
cur mostly by displacing domain walls. Domain bridging 
can be visualized by subtracting two subsequent images 
as shown in Figj2^c). There we see four bridging events 
having roughly the same area, together with a swarm 
of smaller avalanches produced by the rearrangement of 
the surrounding domains. This avalanche triggering phe- 
nomenon is always present when the domain density is 
large. 



B. Characteristic lengthscales 

Our simulations suggest the existence of an upper 
bound for avalanche sizes which appear to be limited 
by the width id of the labyrinthine domains. We es- 
timate the domain width by minimizing Eq.Q with re- 



spect to id assuming a simple form for the magnetization 
field. In particular, we consider stripe domains described 
by m(x,y,£d) = sm(7rx/id) and obtain id = a/j for 
H r = H e = and neglecting the double well potential 
V(m). Figj2|d) compares the domain width estimated 
from simulations at h e = with the theoretical estimate: 
the agreement implies that neglecting the anisotropy in 
the calculation appears to be a reasonable approxima- 
tion. A better estimation has been obtained in Ref. [16j 

A second important length-scale is given by the do- 
main wall width £ w which should induce a lower cutoff 
to the avalanche distribution. If avalanche scaling is due 
to the dynamics of domain walls in a disordered medium 
then for lengthscales below i w the scaling should break 
down since the very concept of domain wall is lost. Below 
this scale the inner structure of the domain wall and the 
rotation of the spins should become important for the 
avalanches, but these features can not be described by 
a simple scalar model. We estimate i w by energy min- 
imization imposing a single domain wall at x = 0. In 
the absence of external field and neglecting dipolar in- 
teraction, Eq.Q is minimized by a magnetization profile 
m(x,y,i w ) — tanh.(x/£ w ) with i w = \f2i xl . The domain 
wall width obtained in simulations is reported in Fig|2jd) 
and compared with the theoretical estimate showing a 
good agreement at least for large 7. 

The third relevant lengthscale in the problem is the 
diameter of the nucleation droplet i n . Bubble domains 
are stabilized by the dipolar i nterac tions and i n has 
been estimated in ideal condition J 18 l 19 l. The theory, how- 
ever, does not include disorder, leading to a overestimate 
of i n . In fact, without disorder (i.e. for D = 0) in 
the parameter range explored here, nucleation is com- 
pletely suppressed and magnetization reversal occurs co- 
herently (i.e. uniformly in the entire sample, without 
domain nucleation/ 2 -^. Hence, we conclude that the nu- 
cleation we observe is induced by disorder, but we were 
not able to estimate i n theoretically. In Fig{2|d) we re- 
port the results of numerical simulations showing that 
i n is slightly larger than id and grows in a similar way 
with a/7. Our theoretical estimates for the character- 
istic lengthscales compare nicely with the experimental 
results of Im et al. on CoCrPt films^ and of Schwarz et 
al. on LaSrMn0 3 nlmPEH. Using K u = 2 x 10 5 J/ra 3 
and A = 3.8 x 10 _12 J/m for a d = 50nm thickness 
CoCrPt film™, we estimate with our model i w = 12nra, 
id = 28nra and i n = 36nra while experiments yield 
i w = 15nm, id = 30nm and i n = 40 nm. For a d = lOOnm 
thickness LaSrMnOs film^, we use K u = 2x 10 4 J/ra 3 and 
A = 1.7 x 10" 12 J/ra and obtain i w 26nm, i d = 63nm 
and i n — 80nm while in experiments we have id = 79nm 
and i n = 73nm, again a very good agreement. 



C. Avalanche statistics 

The statistical analysis of the avalanche distribution is 
performed exactly as in experiments^^. We resolve sin- 
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gle avalanches by computing the difference between con- 
secutive images of the magnetization and identifying the 
avalanche size S as the area of a connected cluster of re- 
versed magnetization, defined by an appropriate thresh- 
old on the phase field (see also supplemental movies 5 
and 6). The probability distributions are obtained by 
logarithmic binning of the measured avalanche sizes, av- 
eraging over several realizations of the disorder field. All 
the measured distributions display a power law behav- 
ior P(S) oc S~ T for at least one decade with the scal- 
ing regime limited at small and large sizes (see Fig{3]). 
Following Ref. 9, we first compute the avalanche distri- 
bution for different values of the applied field along the 
hysteresis loop as shown in Figj3|a) for 7 = 0.7. The 
four displayed distributions span a region ranging from 
the high fields where nucleation takes place to the cen- 
ter of the loop where bridging events appear, passing 
by intermediate regions where avalanches are mainly due 
to domain wall motion. We notice that the distribution 
upper cutoff decreases as we approach the center of the 
loop, where the domain density is maximum. In most 
distributions, we observe a characteristic peak at large S 
that is superimposed to the apparent power law decay. 
We suggest that these characteristic sizes are associated 
with nucleation events at small fields and bridging events 
at larger fields, since by removing those events we observe 
a marked decrease in the peaks amplitudes. A power law 
fit done excluding large values of S yields r = 1.0 in the 
very first stage of the nucleation process when few iso- 
lated droplets expand freely, r = 1.3 when the droplets 
increase their density, r = 1.5 in the intermediate region 
and r = 1.8 in the loop center. While we notice that the 
presence of a peak could strongly bias the estimate of the 
exponents, the values we get are in good agreement with 
the experimental results of Ref. 9 which report a similar 
trend in the variation of r along the loop. 

Next, we distinguish between newly nucleated droplets 
and the expansion of existing domains and measure the 
distributions of nucleated droplet sizes which follows a 
Gaussian distribution (see Fig{3^b)) with a characteristic 
size decreasing with 7. The areas of annihilated domains 
are also Gaussian distributed but their typical sizes is 
independent of a/7 and smaller than those of nucleated 
droplets. This is because, contrary to nucleation, the 
annhilation process does not involve an energy barrier. 
The domain size shrinks continuously until it becomes 
comparable to the domain wall width, when domains be- 
come unstable and disappear. Similar observations was 
made experimentally in Ref. I10|lll In Figj3] (c) we show 
the avalanche size distribution measured in the same cen- 
tral region of the hysteresis loop for different values of 7. 
The avalanche cutoff increases when 7 is decreased, fol- 
lowing the corresponding increase of the domain width 
and confirming that a/7 is the relevant parameter con- 
trolling the size of the scaling regime, depicted by a shad- 
owed region in FigJ5](d). 

Finally Figj3] (d) shows the effect of disorder strength 
on the avalanche size distribution in the central part of 



the loop. The upper cutoff is found to increase with in- 
creasing disorder strength D. A sim ilar result is found 
in domain wall depinning models^ 2 ^, confirming that in 
this regime avalanches are due to domain wall motion. 
The inset of Figj3^d) refers to 7 = 0.7 and shows that, 
even if the domains shape is slightly affected by the dis- 
order strength, is almost independent of D, while the 
nucleation diameter £ n decreases with increasing D. This 
behavior has a simple explanation: to stronger disorder 
corresponds a larger value of the nucleation field H ei but 
higher values of H e lead to a smaller nucleated droplet 
diameter £4. 



IV. DISCUSSION 

In this paper, we have described the process of 
labyrinthine domain formation using a disordered phase- 
field model of magnetic reversal in films with out-of- 
plane anisotropy. The model allows to overcome the 
limitations of spin models, where lattice effects can ef- 
fectively pin domain walls 12 , and at the same time pro- 
vides an alternative to computationally more intense mi- 
cromagnetic equations. We have modeled disorder using 
random-fields in analogy with other well studied hystere- 
sis models 22 , but our main predictions should not change 
substantially with other forms of disorder. To check this, 
we have implemented also random anisotropics 13 find- 
ing qualitatively similar results. The model successfully 
reproduces the main features observed experimentally in 
out-of-plane magnetic thin films^UM: nucleation of nearly 
circular domains, domain growth and branching instabili- 
ties, labyrinthine domains expanding through avalanches 
and finally domain annihilation. Using simple arguments 
we estimate the main characteristic lengthscales of the 
problem finding a good quantitative agreement with ex- 
periments. The domain wall and domain widths, i w 
and id, are shown to delimit the scaling regime for the 
avalanches: the power law distribution of avalanches is 
only observed between lower and upper cutoffs, deter- 
mined by t w and respectively. Such a limited scal- 
ing regime is expected to be the cause for the variabil- 
ity of the exponents reported in the literature for this 
class of materials. The competition between dipolar and 
exchange interactions yield the observed labyrinthine do- 
main structure and at the same time set strict boundaries 
for the scaling region. 

By measuring the avalanche size distribution along the 
hysteresis loop, we find effective exponents that vary con- 
siderably as a function of the applied field in quantitative 
agreement with experimental results for CoCrPt films^. 
We tend to attribute this variability in the exponents 
to a bias in the fit due to the upper and lower cutoff. 
The most reliable value we obtain is r c± 1.5 observed 
in the intermediate regime where nucleation and bridg- 
ing events are not valid. This value is very close to the 
prediction of mean-field theory that is found to be accu- 
rate for bulk samples 2 -^, although we have no theoretical 
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FIG. 3: Panel (a) shows the avalanche probability distribution for 7 = 0.7 measured in different regions of the hysteresis loop. 
The region is indicated in terms of the applied field h e interval together with the fitted power law exponents r. The curves 
are shifted to help comparison. Logarithmic binning has been adopted. Panel (b) shows the probability distribution of the 
areas of nucleated droplets and for the areas observed just before domain annihilation. The inset shows the size of the very 
first nucleated droplet for different simulations with growing 7. Panel (c) shows the avalanche size distribution for different 7 
measured in the central region of the loop with h e ranging between —0.5 and 0.5. Panel (d) shows the probability distribution 
of avalanche sizes for 7 = 0.5 at different disorder strength D. The inset highlights the change in the domain morphology 
induced by the disorder strengthening. 



argument to justify it for thin films. The values we find 
are, however, in net contrast with the one reported for 
LaSrMnOs film (r ~ 0.5) 10 | 11 [ We notice, however, that 
the result is due to avalanches that occur for lengthscale 
smaller than the domain wall width £ w . For this reason, 
they can not be captured by our model that does not in- 
volve spin rotation and other small scale details. Yet, we 
have reanalyzed the data of Ref. 110)111 and found that 
on lengthscales larger than i w one sees a larger exponent 
although its determination is masked by the large scale 
cutoff. 

The results of our numerical simulations yield some 
predictions that could be tested experimentally. Fig. 
[3^c) shows that the upper cutoff of the avalanche size 
distribution increases with 7. Experimentally it would 
be possible to tune 7 by analyzing magnetic films of 



different thickness. Similarly, Fig. [3jd) reports a similar 
effect due to the strength of the disorder, which could be 
changed experimentally by irradiating samples, as shown 
in Ref 15 . Given the quantitative agreement between 
our model estimates and the experimentally measured 
lengthscales, we are confident that this approach should 
be successful. 
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